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Abstract 

This work concerns the simulation of compressible multi-material fluid 
flows and follows the method FVCF-NIP described in the former paper 
This Cell-centered Finite Volume method is totally Eulerian since the mesh is 
not moving and a sharp interface, separating two materials, evolves through 
the grid. A sliding boundary condition is enforced at the interface and mass, 
momentum and total energy are conserved. Although this former method 
performs well on ID test cases, the interface reconstruction suffers of poor 
accuracy in conserving shapes for instance in linear advection. This situation 
leads to spurious instabilities of the interface. The method Enhanced-NIP 
presented in the present paper cures an inconsistency in the former NIP 
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method that improves strikingly the results. It takes advantage of a more 
consistent description of the interface in the numerical scheme. Results for 
linear advection and compressible Euler equations for inviscid fluids are pre- 
sented to assess the benefits of this new method. 
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1. Introduction 



The two-material compressible hydrodynamics equations (Euler equa- 
tions) are considered in this work. The flow regime is such that molecular 
viscosity within materials is neglected: materials are supposed immiscible 
and separated by sharp interfaces, with perfect sliding between materials. 
Each material is characterized by its own equation of state (EOS). 

The formalism of finite volume methods is close to the mechanical view- 
point, and generic for different types of physical models. Thus, it might be 
easier to add such models; surface tension or turbulent diffusion for instance. 
The discretization order is limited, but this method is accurate to simulate 
hydrodynamic shock waves, because of the consistency between numerical 
treatment and mechanics. 

The extension of Eulerian schemes to multi-material fluid flows can be ob- 
tained by various techniques. One is to introduce the cell mass fraction c a of 
material a and let it evolve according to material velocity. The cell is called 
pure if a material a satisfies c a = 1 and is called mixed if c a g]0, 1[. Pure 
cells filled by material a are calculated in the same manner as for the single 
material method. Mixed cell evolution is computed using a mixing equation 
of state that takes into account material mass fractions, see e.g. l|. One 
drawback of this approach is the numerical diffusion of the interface that ob- 
viates sharp interface capturing. It turns out that for some applications, this 
drawback is not acceptable since the diffusion of one material into another 
one will correspond to a different physics. For example the two material 
could react when a molecular mixture is formed. Hence such a diffusion 
should occur only for physical reasons and not for numerical ones. 
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In the case of sharp interface capturing methods, the interface is approx- 
imated in a mixed cell by a segment by most authors. However more com- 
plex curves than straight line or more complex theory (see 6[ for instance) 
might be used. A famous method using sharp interface reconstruction is 
the L agr ange-f- Remap Finite Volume scheme, initiated in [9] and further im- 
proved in 10j. It belongs to the family of so called Volume of Fluid (VOF) 
methods. The first step of this method is a Lagrangian scheme, resulting 
in a mesh displacement with material velocity. The second step is a multi- 
material remapping of Lagrangian mesh onto the original Eulerian mesh, by 
exchanging volume fluxes between cells related to the Lagrangian motion of 
cell edges. The new interface position in mixed cells is determined using 
the partial volumes of the materials and the interface normal vector. The 
later is calculated using volume fractions from neighboring cells. Thus the 
ratio of each material in volume fluxes is deduced from the multi-material 
remapping. Some methods with the same kind of operator splitting are used 
for incompressible multi-material fluid flows as in |8j. These methods pro- 
vide sharp interface between materials and discontinuous quantities in mixed 
cells, allowing large deformations and transient flows. In this context, the 
drawback of these Lagrange+Remap methods is the limited accuracy of the 
underlying single phase scheme due to diffusion induced by the remapping 
step. Moreover, more complex physics at material interfaces such as sliding 
effects, is not possible. 

The FVCF scheme (Finite Volume with Characteristic Flux) has been 



introduced in 



71 



for simulating single phase compressible flows or multi- 



phase models without sharp interface capturing. The method described in 
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5], so called NIP method (Natural Interface Positioning), is an add-on to 
the FVCF method in order to deal with multi-material fluid flows with sharp 
interface capturing. It is a cell centered totally Eulerian scheme, in which ma- 
terial interfaces are represented by a discontinuous piecewise linear curve. A 
treatment for interface evolution is proposed on Cartesian structured meshes 
which is locally conservative in mass, momentum and total energy and allow 
the materials to slide on each others. Discrete conservation laws are writ- 
ten on partial volumes as well as on pure cells, considering the interface in 
the cell as a moving boundary without any diffusion between materials. A 
specific data structure called condensate is introduced in order to write a 
finite volume scheme even when the considered volume is made of moving 
boundaries, i.e. interfaces. This treatment includes an explicit computation 
of pressure and velocity at interfaces. 

In [5( are shown 2D results illustrating the capability of the method to 
deal with perfect sliding, high pressure ratios and high density ratios. This 
former method however produces non satisfactory results in the context of 
advection of geometrical shapes especially when dealing with low Mach num- 
bers. It is however a classical misbehavior of most of advection and recon- 
struction methods which have a tendency to destroy the shape of advected 
objects due to numerical approximations. However, this former method gives 
very poor results when advecting geometrical shapes especially when dealing 
with low Mach number flows. In this work we propose a new method called 
ENIP (Enhanced NIP) that is an improvement of the NIP method by a more 
accurate treatment of condensates. On a very simple example: the advection 
of a square, an inconsistency in the NIP interface reconstruction method will 



5 



be exhibited. We will then introduce ENIP that cures this situation. Nu- 
merical examples are presented in the last Section to assess the validity and 
efficiency of this new approach. 



2. FVCF-ENIP: Finite Volume Characteristics Flux with Enhanced 
Natural Interface Positioning technique 

2.1. Governing equations 

The model addressed in this work is the compressible Euler equations in 
space dimension d that can be written in a conservative form as follows: 

^( p ) + div(H = o, (i) 

d 

— (pu) + div(pu®u + pl) = 0, (2) 
ot 

^(pE) + div {{pE + p)u) = 0, (3) 

where p denotes the density, u £ R d the velocity field, p the pressure, 
E = e + \u\ 2 /2 the specific total energy and e the specific internal energy. 
An equation of state of the form EOS(p, e,p) = or p = p(p, e) is provided 
in order to close the system. 

Let us consider a generic conservative form with V = (p, pu, pEf the un- 
known vector of conservative variables and flux F is a matrix valued function 
defined as: 

F : R d+2 — > R d+2 x R d 

, ^ (4) 

for all direction n £ R d , F(V) ■ n is given in terms of V by: 

F(V) ■ n = (p(u • n) , pu{u • n) + pn, (pE + p) (u ■ n)) . (5) 
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The compressible Euler equations (fUE]) can then be rewritten as: 

d t V + dwF(V) = 0. (6) 

2.2. FVCF: Single material scheme 

FVCF method uses a directional splitting on Cartesian structured meshes. 
The method is thus detailed for only one generic direction denoted by x. In 
d dimensions of space, the algorithm described for direction x has to be 
replicated d times, one for each direction. However, this directional splitting 
does not modify at all the underlying single material scheme FVCF for pure 
cells. In 2D: 

- variables at t n ' x are calculated from those at t n by the x direction step, 

- variables at t n+1 are calculated from those at t n ' x by the y direction 
step. 

v n,x _ y n 

Vol, * — 1 + A. + = 0, (7) 

Vol, 1 At 1 +A y (r d + 4> n u) = 0, (8) 

where the cell volume is Vol,, the cell face area are A x and A y respectively 
normal to x and y directions, up, down, right and left direction fluxes (/>", 
0^, 0", 0™ calculated with respect of the outgoing normal direction of 
cell face Td in direction d using variables at time t n , i.e. 

^ = j-J F ( yn ) ■ n * ds - ( 9 ) 

This flux is 
scribed in 



urther approximated using the finite volume scheme FVCF de- 



2.3. FVCF-NIP: Multi-material scheme 

One considers multi-material flows. The subcell model addressed here for 
the multi-material representation is a cell C of volume Vole* containing n m 
different materials, each of them filling a partial volume Vol^ such that 

' ' ■in 

^Vol£ = Vol c . (10) 

fe=i 

Cell C is referred to as pure if n m = 1, and as mixed if n m > 1. The interfaces 
in mixed cells are approximated by segments separating materials into two 
partial volumes which are pure on both sides of the interface. 
A partial volume cell-centered variable vector Vk = (pk, PkUk, PkEkY and an 
equation of state EOSk(pk, £k,Pk) = are also associated with each material 
labeled by k < n m in the mixed cell. 

FVCF-NIP method uses a directional splitting scheme for the interface evo- 
lution without loosing the accuracy of the Eulerian scheme in the bulk of 
materials. Consequently this scheme is restricted to structured Cartesian 
mesh. 

The multi-material extension proposed in [5( considers the finite volume 
scheme (EHHJ) on each partial volume in a mixed cell. The obtained scheme 
is conservative by construction and is constrained with the same CFL con- 
dition as the single material schemqj. NIP method consists in removing cell 
edges when this cell contains an interface. Therefore each partial volume is 
merged with the neighbor pure cells filled with the same material, see Figure 
[TJ Variables in these enlarged partial volumes are obtained by writing the 



1 Without such a special treatment the time step would be constrained by the smallest 
partial volume, which is arbitrarily small. 
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f" Pure cell Pure cell Mixed cell Pure cell Pure cell 




Interlace 



f™* ■* Condensate 




t Pure cell Pure cell Pure cell Mixed cell Pure cell 




Interlace 



Figure 1: Sketch of a condensate. Evolution of an interface through a cell edge during one 
time step. Mixed and pure neighbor cells are merged to obtain the so called condensate 
at fictitious time t n *. Interface evolution is performed within this condensate from t n * to 
t n+1 * This condensate is then split back into Eulerian cells. 



conservation laws on the merged volumes 



Vol! = Voli+Vol^i, (11) 



Vol 2 = Vol 2 +Vol pure2 , (12) 
then on the conserved variables 

y _ Volj Vj + Vblpure 1 Vpure 1 

1 voir 

y Vol 2 V 2 + Volpure 2 Vpure 2 

2 VoT 2 

This set of cells is associated with its left and right single material fluxes (fie 
and <fi r . Internal cell edges are forgotten, considering only enlarged volumes 



Voli and Vol 2 and averaged variables V\ and V 2 , separated by an interface; 
this system is called a condensate. 

Actually, this numerical strategy consists in condense neighboring mixed 
cells in one direction of the Cartesian mesh, in which interfaces are con- 
sidered as mono dimensional objects, namely they are considered vertical 
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during x direction step and horizontal during y direction step. A condensate 
then contains layers of successive different materials that are separated by 
straight interfaces. The thickness of these layers is calculated through vol- 
ume conservation. The ordering of layers is given by the 2D description from 
the previous time step. It is determined thanks to the volume fractions of 
neighboring cells. The layer evolution is calculated in a Lagrangian fashion 
which implies that layers can be as thin as partial volumes are small. Once 
quantities and interface positions inside the condensate are known at time 
t n+ , they are remapped back onto the original Eulerian mesh. Finally a 2D 



normal in each mixed cell is computed as described in 10]: the method is 
based on an approximation of the gradient of the volume fraction function 
in mixed cells. It provides the normal to materials interface in each cell 
that is further used to locate materials within mixed cells. The numerical 
scheme used in a condensate is presented in great details in [5] and we omit 
this description in this work and rather focus on the interface reconstruction 
method. 



As shown in 



5] this numerical method has several attractive properties 



as conservation and perfect sliding of materials as instance. Moreover At is 
not restricted by small partial volume thanks to a tight control of density 
and pressure The numerical experiments carried out in 7, 2, 5, J] have 
confirmed the efficiency of such a method for compressible multi-material 
computation. Although very promising, the method suffers from the way 
interfaces are dealt with. 

In order to illustrate the interface reconstruction method NIP let us con- 
sider a square like interface cutting the Eulerian cells, as in Figure 121(A) . 
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These interfaces are indeed defined by their normals within each cell. NIP 
method consists of the following steps assuming the condensate is in the x 
direction: 

• Representation Figure E|-(B). The representation step can be seen as 
the way of determining on which side (left or right) of the mixed cell 
the material is to be put. This is done by comparing the direction of 
the interface normal at time t n with the vertical direction. 

• Condensate construction Figure [2]- (C). The construction of the conden- 
sate consists in discarding any cell edges in the mixed cells considered. 
Then the partial volumes of the same contiguous materials are glued 
together into so called condensate layers. As instance cell 2 and 3 dark 
materials are merged into one stand-alone layer with associated volume 
averaged values. 

• Condensate evolution Figure [2]-(D). The condensate layers evolution 
is computed from t n to t n+1 thanks to the numerical scheme devel- 
oped in [5(. In short, each vertical interface is assigned a velocity and, 
consequently, a new position of each layer within the condensate is de- 
termined in a Lagrangian way. Any conserved variable is computed 
accordingly. 

• Reconstruction Figure [2]-(E). This phase consists in "guessing" the 
shape of each layer in the condensate before remapping. The recon- 
struction phase was not originally considered as a true phase of the 
algorithm as the author used the same shapes as the ones produced in 
phase Condensate construction, i.e. only vertical interfaces. 

11 
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Figure 2: NIP method — (A) Situation at t n with real materials geometry, interfaces 
and normals to them. (B) Representation of partial volumes at t n . (C) Construction of 
a condensate at t n by merging layers of contiguous partial volumes of the same material. 
(D) Evolution of condensate in a Lagrangian fashion during At. (E) Condensate re- 
construction at t n+1 . (F) Condensate projection/remapping from Lagrangian mesh onto 
original mesh. 
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• Projection Figure |2]-(F). The projection step consists in remapping the 
shapes obtained from the reconstruction phase onto the Eulerian grid. 
This step produces updated partial volumes in mixed cells. Volume 
fractions are deduced. 

When all mixed cells in the domain are treated for direction x, the interface 
normals are computed using the updated volume fractions. This concludes 
the system evolution in direction x, as we are back to a similar situation as 
the one described in Figure [2]- (A). 

In the case where the normal is almost vertical, positioning the material on 
either side of the cell might be, at least inaccurate, or, worse, incorrect. Fur- 
thermore the reconstruction phase is here clearly inconsistent: the interfaces 
are initially horizontal in cell 3 and 4 (Figure [2]- (A)), while in the Recon- 
struction Figure |2]-(E) and in the Projection Figure IH-(F) phase interfaces 
are set vertical for any initial geometry. This situation of a horizontal inter- 
face is the worst case, but it illustrates the lack of geometrical consistency of 
NIP. This inaccurate reconstruction step leads to a lack of accuracy of the 
volume fractions obtained after the remapping step. Ultimately, it impacts 
the whole numerical method in any advection process. 

As an illustration let us consider the diagonal advection of a square back and 
forth as shown in Figure [3J We omit the exhaustive description of this test 
as it will be done in the numerical Section of this paper. On the right panel 
it is obvious that the shape of the square is not well approximated. More 
important the horizontal and vertical edges of the square do not remain so. 
This behaviour is less pronounced if one refines the mesh but still remains. 
Our goal is to improve the reconstruction step so that the new method, de- 
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Figure 3: Diagonal advection of a square by NIP method — Left: Exact configuration — 
Right: Numerical configuration obtained by NIP. 



noted ENIP standing for Enhanced Natural Interface Positioning, cure this 
geometrical inconsistency during the advection phase of the algorithm. 

24. FVCF-ENIP 

The main idea of the new interface reconstruction method ENIP emanates 
from the following remarks: 

1. At time t n any interface normal in mixed cell i denoted nj is known. It 
is used to locate the partial volumes within cell i when the condensate 
is constructed (phase (B) and (C) of Figure [2]). However is never 
taken into account in the reconstruction and projection phases (E) and 
(F) from the same figure. 

2. Any layer of the condensate evolves as a Lagrangian object in the orig- 
inal method. Consequently the cell faces could evolve in an almost 
Lagrangian manner within this condensate. This makes possible to 
conserve the initial geometry of partial volumes during this Lagrangian 
motion. 
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Therefore ENIP modifies several steps of NIP as depicted in Figure HI Once 
a patch of neighbor mixed cells in x direction^ are agglomerated, The same 
five steps as for NIP method are performed. The first two steps are kept 
unmodified. The last three are modified as described in the following. 

2.4-1- Lagrangian Condensate evolution step 

Cell interface Lagrangian velocity. After the condensate at t n is constructed, 
each layer labeled c is located thanks to the left and right interface position 
respectively called x~,x^. The numerical scheme provides the layer evolu- 
tion, and as a by-product, the velocity of these interface positions, u~,u^ 
are given by 

x c --"+ 1 = x ~ + At u~, xp n+1 = x+ + At u+. (15) 

We make the following fundamental linear displacement assumption: The 
velocity linearly varies within any layer, see Figure [5] for a sketch. This 
assumption implies that any point +] characterized by its ID 

barycentric coordinates 

ry> ~ \~ ry . ry , ry 

\ = —r r> \ = —r r> (16) 

ry\^ ry | ry 

moves to location 



X ; 



" +J - K^ +1 + \ + *+' n+1 = Xi + At (X[u~ + Xfut) • (17) 



Then the point velocity is naturally set to Ui = X[u~ + \~l u t ■ Using this 
previous formula one can associate a "Lagrangian" velocity to any cell inter- 
face. As instance in Figure ffl(C) cell interface located at x" moves to the 



2 The y direction is treated likewise. 
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Figure 4: ENIP method — (A) Situation at t n with real material, interfaces and nor- 
mals to them. (B) Representation of material at t n . (C) Construction of a condensate 
at t n by merging of mixed cells leading to layers of contiguous pieces of the same mate- 
rial. (D) Evolution of condensate in a Lagrangian fashion during At. Determine layer 
compression rates <5Vol^ through the evolution of Lagrangian cells during At. (E) Con- 
densate reconstruction at t n+1 using interface normals defined at t n . (F) Condensate 
projection/remapping from Lagrangian mesh onto Eulerian mesh. 
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Figure 5: Sketch of linear displacement assumption — Displacement velocity varies lin- 
early between the layer interface velocities (x in color) computed by the numerical scheme. 
The cell interface velocity (•) is interpolated. The two top rows represent the evolution of 
a condensate in the x direction from t n to t n+1 . 



position x™ +1 = x™ + At Ui with Ui being the linear combination between 
u~ and m+ via the barycentric coordinates of point Xi in [x~; With the 
same formula one gets x™+{ = xf +1 + A u~ +1 in the next layer as = x~ +l . 

Compression/ expansion rates. The global rate of compression/expansion in 
layer c during At is given by 

ry~\~ i^ - ) - 1 ry " j^H - ! n i ~T~ q i 

5V ol c = c - c = 1 + At 1 c . (18) 

iy~j ry ry~\~ iy 

The linearity assumption provides a simple way to determine the rates of 
compression/expansion at left/right of a point 

T n+1 _ ™-,n+l +,n+l _ T n+1 

5Vol: = Xi - Xc , 5Vol+ = ^ (19) 

that fulfil cWbl~ + 5Vol^ = <5Vol c . Moreover the substitution of x™ +1 in the 
previous equations yields 

5V oL7 = X :~ X < + At U :~ U < = X+ + At Ul ~ U ~ c , (20) 

rf~ir />"■ iy~\~ rn iy ~T~ ry — 

c c c c c c 
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where Ui—u c = (\ u c +^t u t)~ u c = ^t( u c~ u c )> therefore the compression 
rates simply writes 

5Vol c ~ = A+ ( 1 + At U {~ Uc ) = A+5Vol c , (21) 

SVolt = K (l + A ^+I^_ ) = \ r ^Vol c . (22) 

Each SVol^ or 5Vol~ is associated to a unique Eulerian cell; as instance in 
Figure IU 5 Vol" is associated to cell 2, <5Vol,t~ to cell 3, ^VoLt^ to cell 4 
and so on. Therefore <5Vol^ provides de facto the compression/expansion 
of the partial volume originating from its associated Eulerian cell motion. 
Furthermore, as any Eulerian mixed cell i possesses a unique normal denoted 
nj, this last is associated to the corresponding partial volume <5VoFr; this 
normal is consequently labeled nf. These rates are then used to reconstruct 
the material topology into the Lagrangian cell. 

2.4-2. Reconstruction step 

The Lagrangian cell i + 1/2 at t n+l the interfaces of which moved as 



x 



Xi + At m, x™^ 1 = x i+ i + At u i+1 , (23) 



changed its volume as 

SVok+1/2 = pi = *_SL^?1 = , + toOaz*. (24) 

The velocity Ui depends on u~_ 1: u^_ x and Ui+\ depends on u~ , u+. Moreover 
u*_ x = u~ by definition. 

The second fundamental assumption states that the interface normals do 
not change their direction during their Lagrangian evolution. The goal is to 
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locate the partial volume into the Lagrangian cell at t n+1 and construct the 
linear interface, knowing its normal nf. Necessarily this partial volume is 
either in contact with cell interface Xi (superscript +) or Xi + \ (superscript 
— ). Its volume at t n+1 is given by 

y±,n + l = y± JVol T = y± + ^ Xf{u + _ (35) 



If V^ n+1 < V^ 2 then there exists a unique line oriented by the normal nf 
and separating the cell volume into two sub-volumes V^' n+1 and (V^y 2 — 
y±,n+i^ reS p ec tively by the PLIC ("Piecewise Linear Interface Construction" 



10l |) method. As the displacement velocity u(x) is supposed to be piecewise 
linear (by the first assumption see Figure [5]), then, if Xi < x~ < Xi + \ one 
deduces x™ +1 < x~' n+1 < x^{. Therefore the sub- volume at t n+1 is strictly 
included into the Lagrangian cell volume V^L. This phase is depicted in 
Figure ll-(E) 

2.4-3. Projection step 

The projection step performs the exact intersection between the La- 
grangian condensate obtained after the reconstruction step in Figure H]-(E) 
and the Eulerian mesh (bold line squares in Figure 0]-( A)). This step is de- 
picted in Figure Hl-(F). The exact intersection consists in projecting each 
partial volume that is accurately located into the condensate, onto some Eu- 
lerian fixed cell(s). As instance in Figure IH-(F) the first partial volume is 
projected onto Eulerian cells 2 (green cell) and 3 (red cell). Contrarily the 
last partial volume is totally projected into Eulerian cell 5 (brown cell). This 
projection provides the quantity of material per Eulerian cell, or, equivalently 
its volume fraction. 
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Once volume fractions in the mixed cells are updated though the evolution 
of condensates, 2D normals are computed using the same technique as in 
original NIP method. 



3. Numerical results 

In this Section we present a set of test cases to assess the efficiency of 
the approach described in the previous Sections. First, one validates the 
technique on pure advection test cases that often present excessive smearing 
of interfaces due to the numerical inaccuracy embedded into the scheme. 
A square shaped object is advected with constant velocity in a diagonal 
direction in a first test, then into a rotating flow. Finally an hydrodynamics 
test case is presented. 

3.1. Advection context 

An initial square [0.1; 0.1] x [0.2; 0.2] is located into the domain VL = [0 : 
0.4] x [0; 0.6]. The density into the square is set to po(x) = 1 whereas it is set 
to po(x) = outside. In the pure advection context this square shape should 
be perfectly conserved through the equation 

l p + u l p + v T/ = a - (26) 

where (u, v) is a constant velocity field. The exact solution at any point x 
and any time t is p ex (x, y,t) — po(x — u t,y — v t). If the numerical method 
provides an approximated solution called pf in cell % at time t n then the error 
in L a norm is evaluated by (a = 1, 2) 

EM - p ex (xi,t n )\ a 



Y,i\p ex {xi,-t n )\ a 
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(27) 




* 



» ■ 



Figure 6: Advection of a square (zoom around the exact position of the initial and final 
square) — Top-left: exact solution — Top-right: classical NIP with a 60 x 60 mesh - 
Bottom-left: classical NIP with a 120 x 120 mesh — Bottom-right: ENIP with a 60 x 60 
mesh. 



The first test consists in advecting the square with the constant velocity 
field u — 1, v — 3 up to the time t = 0.1 then reversing the advection 
field by setting u = —1, v = —3 up to final time t = 0.2 so that the final 
configuration perfectly fits the initial one. Any method (NIP and ENIP 
included) introduces some error that we intend to measure with this test. 
In Figure |6] are shown the exact solution (top-left) and the results obtained 
with a 60 x 60 mesh for NIP (top-right) and ENIP (bottom-right). ENIP is 
visibly able to preserve the shape of the square whereas NIP is not. A mesh 
refinement of NIP computation (120 x 120 mesh for the bottom-left panel) 
does not improve the situation. In tabled] we gather the errors for the L 1; 
L 2 norms for successively refined meshes for the NIP and the proposed ENIP 



21 



Ax = Ay 


Lt NIP 


L x ENIP 


L 2 NIP 


L 2 ENIP 


0.02 
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0.0133 
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0.01 


0.339 


0.111 


0.284 


0.053 


0.005 


0.221 


0.042 


0.195 


0.017 


0.0033 


0.155 


0.025 


0.138 


0.010 



Table 1: Error in L\, L2 norms for the advection problem — NIP versus ENIP methods. 
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Figure 7: Convergence of ENIP vs NIP for a pure advection problem. The log of the L2 
error is displayed as a function of the log of Ax. 

method on this advection problem. Systematically ENIP over-tops NIP. In 
Figure [7] we display the log- log scale results for the error in L 2 norm for both 
methods showing the improvement gained by ENIP; indeed the slope which 
represents a measure of the numerical order of convergence is improved by a 
factor 2.5 (0.6 for NIP and 1.5 for ENIP). 

The next test consists in the rigid rotation of a square [0.06; 0.46] x 
[0.3; 0.7] (density 1) into the unit square domain, see Figure M top- left panel. 
A 100 x 100 uniform mesh is considered and the rotation is given by the 
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Figure 8: Top-left: Sketch of the rigid rotation of a square — Top-right: After 5/8 of 
the full rotation — Bottom-left: After one full rotation — Bottom-right: After three full 
rotations. 

velocity field 

u = -100(y-0.5), v = 100(x-0.5). 

In Figure M we display the density after 5/8 of the full rotation, after one 
and three rotations. The square shape is almost preserved. Contrarily the 
classical NIP method would totally lose the shape after one rotation. 

3.2. Hydrodynamics context 

We run an idealized 2D test case that corresponds to the free drop of 

n 

a liquid rectangle within a 2D rectangular tank filled with gas |4|. This 
context is inspired by the problem of sloshing that may appear in the tanks 
of Liquid-Natural-Gas (LNG) carriers. The study focuses on the ability for 
the numerical simulations to take properly into account the physics that is 
of major importance during the liquid impact such as the escape of the gas 
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underneath and its compression. As a strong sliding process occurs between 
the compressed gas and the falling liquid. The ability of the method to 
properly deal with sliding conditions at the interface has a major effect on 
the final numerical compression and shape of the trapped air. This has 
ultimately a strong influence on the impact pressure. 

The test case consists in a domain Q = [0.0; 0.0] x [10m; 15m] filled with air. 
The liquid is initially at rest in the rectangle [0; 2] x [5; 10] and is falling under 
the gravity that is pointing downward with magnitude g = 9.81m. s~ 2 . A 
free fall of the liquid into vacuum would impact at ^i m pact = 0-64s however 
due to the presence of the gas this theoretical value is not correct for our 
simulation however some critical phenomena still occur in the vicinity of 
this time. As instance around ^i m p ac t a pocket of gas is trapped under 
the falling liquid and this strongly impacts the numerical impact pressure 
by decelerating and damping the free fall of the liquid. Therefore a good 
interface reconstruction method should qualitatively improve the numerical 
results. One considers a mesh made of 100 x 150 uniform cells on the domain. 
One shows the results for NIP and ENIP at time t = 0.6s Figure |9]-(a)-(b) 
and t = 0.64s in Figure[9]-(c)-(d). The classical NIP method was already able 
to deal with such sliding effects. However the interface reconstruction method 
employed is not accurate and stable enough to be free of oscillation that 
one suspects to be only a numerical artifacts (see panels (a-c)). Contrarily 
the new reconstruction method ENIP on this very same test case is able to 
produce a smooth interface that permits to obtain a more realistic simulation. 
Indeed this simulation prominently displays the fact that the "bubbling" 
effects of NIP is of pure numerical origin and that ENIP cures this drawback. 
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(c) 



(d) 



Figure 9: — (a) NIP results at t = 0.6 (full view and zoom on the impact zone) — (b) 
ENIP results at t = 0.6 — (c) NIP results at t = 0.64 — (d) ENIP results at t = 0.64 . 



4. Conclusion and perspectives 

This paper deals with the improvement of the so-called NIP (Natural In- 
terface Positioning) method. The NIP method described in js] is an add-on 
to the FVCF method in order to treat multi-material fluid flows uses the 
concept of condensate. A condensate is the association of contiguous mixed 
cells in either x or y direction. They are further treated as an entity to 
make possible the treatment of each mixed cell taken individually. NIP is 
the method based on the following steps: Representation, Condensate con- 
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struction, Condensate evolution, Reconstruction, and Projection. The present 
paper points the weakness of the NIP method in pure advection context and, 
consequently, in a full multi-material hydrodynamics one. An enhanced NIP 
method is proposed (ENIP). It modifies several of the previous listed steps. 
More precisely the condensate is assumed to evolved in an almost-Lagrangian 
fashion. The reconstruction step assumes that the condensate keeps the same 
form modulo some expansion/compression that the numerical scheme already 
provides. So the displacement of the condensate is performed either with the 
true computed velocity or with an interpolation of it. In fine the condensate 
preserves its topology contrarily to the original NIP method for which the 
condensate has no recollection of its shape from the beginning of the time 
step. 

The capability of the full numerical method is now dramatically improved 
as seen on advection test cases (advection and rigid rotation of a square). 
Moreover we ran ENIP on a difficult mutli-material hydrodynamics tests 
simulating the free drop of a liquid rectangle within a 2D rectangular tank 
filled with gas in the context of sloshing that may appear in the tanks of 
Liquid- Natural- Gas carrier (see [J]). The accuracy, stability and robustness 
of the ENIP method is clearly seen especially at the time some air is trapped 
under the water. In the near future we plan to investigate the evolution of 
this method to the case of mixed cells with more than two materials. In this 
case the only difficulty lays in the positioning of the different materials in 
the cell, but their evolution within the condensate follows exactly the same 
algorithm ENIP with no modification of the numerical scheme. We also plan 
to investigate the evolution of the method in 3D. 
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